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Abstract 

The least-squares finite element method (LSFEM) is applied to electro- 
magnetic scattering and radar cross section (RCS) calculations. In contrast 
to most existing numerical approaches, in which divergence-free constraints 
are omitted, the LSFEM directly incorporates two divergence equations in 
the discretization process. The importance of including the divergence equa- 
tions is demonstrated by showing that otherwise spurious solutions with large 
divergence occur near the scatterers. The LSFEM is based on unstructured 
grids and possesses full flexibility in handling complex geometry and local re- 
finement. Moreover, the LSFEM does not require any special handling, such 
as upwinding, staggered grids, artificial dissipation, flux-differencing, etc. Im- 
plicit time discretization is used and the scheme is unconditionally stable. By 
using a matrix-free iterative method, the computational cost and memory re- 
quirement for the present scheme is competitive with other approaches. The 
accuracy of the LSFEM is verified by several benchmark test problems. 


1 Introduction 

The most widely used numerical method for solving the time-dependent Maxwell 
equations has been the finite-difference time-domain (FDTD) scheme developed by 
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Yee [44], and extensively utilized and refined by Taflove et al. [38] and Kunz and 
Luebbers [21], as well as others (see a recent survey by Shlager and Schneider [35]). 
In the FDTD, only two Maxwell curl equations are solved by using the explicit 
time-marching scheme and the central difference approximation. According to Yee’s 
scheme, the electric and magnetic fields axe located at different nodes and computed 
alternatively in time. In other words, the computational mesh for the FDTD has 
to be structured, and staggered in both time and space to prevent an oscillatory 
solution from developing as the calculations proceed. 

The numerical techniques developed in computational fluid dynamics, based on the 
conservation laws, have also been applied to solve the time-dependent Maxwell curl 
equations, including finite-difference (FDTD), finite volume (FVTD) and finite ele- 
ment methods (FEM). All these approaches require some sophisticated treatment, 
such as upwinding with characteristic-based flux-differencing by Shankar et al. [34] 
and Shang [31], staggered grids with artificial dissipation by Noack and Anderson 
[28], or the Taylor- Galerkin method with flux corrected transport by Ambrosiano et 
al. [2], It should also be noted that, only for high-speed compressible flow problems, 
in which shocks occur due to the nonlinearity of the aerodynamics equations, the 
conservation form is more appropriate, since this formulation often provides a better 
prediction of the location and strength of the shocks. However, the Maxwell equa- 
tions in electromagnetics are linear. Except at material interfaces, the solution for 
such a system of equations is C° continuous, i.e., it does not contain sharp, shock-like 
discontinuities. Therefore, the numerical simulation of problems in computational 
electromagnetics should not have to resort to the conservation form of equations 
and the above-mentioned complicated techniques which are closely associated with 
shock-capturing. 

The FEM has become the dominate method for numerical solution of problems in 
static electromagnetics, eddy currents and waveguides, see, e.g., Silvester and Pelosi 
[36]. However, application of the node-based FEM to time dependent microwave 
analysis is rare. Mei et al. [7], Madsen and Ziolkowski [23], Choi et al. [10], Wong 
et al. [42] and Morgan et al. [25] are among the few papers which can be found in 
this regard. 

In all the above mentioned time-domain approaches, only two curl equations are 
solved, and the divergence equations are neglected. For frequency-domain computa- 
tion, such as numerical simulation of waveguide and eddy current problems, it is well 
known that if the divergence-free constraint is neither forced nor implied, spurious 
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modes occur and they give particularly high values of V • B, see e.g., Rahman et al. 
[30], Davies [12], Jin [19] and references therein. However, very few have realized 
that ignoring the divergence equations in the time-domain computation also leads 
to spurious solutions. Wu and Jiang [43] first reported some evidence that clearly 
shows the significant violation of the divergence-free condition near the boundary of 
scatters in the solutions of the curl equations only. Kangro and Nicolaides [27] also 
found that time marching solutions of the Maxwell curl equations are contaminated 
by spurious stationary components and gave a simple technique to remove them. 
Indeed, it is commonly believed that the divergence-free constrains in the Maxwell 
equations are redundant, see for example, Kunz and Luebbers ([21], pg. 11). This is 
usually attributed to the following consideration: taking the divergence of the Fara- 
day and Ampere laws, one finds that these divergence-free conditions are satisfied for 
all time if they are satisfied initially. However, it is not easy to satisfy them initially 
by use of the conventional methods, since the inclusion of them will make the system 
of partial differential equations “overspecified” or “overdetermined”, i.e., there will 
be more equations than unknowns. It is common practice to assume that the initial 
field intensities are zero throughout the domain, and that the boundary conditions 
on the surface of the scatterer are correctly given. In this case, the divergence-free 
condition is significantly violated near the scatterer surface in the first time step 
of computation. Therefore, there is no guarantee that the computed fields will be 
divergence-free. However, in many scattered field calculations no attention was di- 
rected to the verification of whether the computed field intensities actually satisfied 
divergence-free conditions. A few exceptions to this are the work by Shang and 
Gaitonde [32] and by Ambrosiano, et al. [2], in which the values of the divergence 
of the computed field were numerically examined. It is not difficult to prove that 
divergence-free conditions are often violated. For example, for two-dimensional TM 
wave problems, one just needs to show the scattered H x and H y contours. Unfor- 
tunately, such results are rarely found in computational electromagnetics literature. 
An exception is the paper by Vinh et al. [40]. In Figs. 2-b and 2-c of that paper, 
a non-physical solution near the circular cylinder can clearly be observed. Such 
spurious solutions are a direct result of neglecting the divergence equations, as will 
be further demonstrated in Section 4.4 of the current report. 

It should be mentioned that Assous et al. [3], followed by Sonnendriicker et al. [37], 
frilly realized the importance of the divergence equations in the time-domain Maxwell 
equations and correctly worked with second-order wave equations. For a special case, 
they supplemented two divergence equations and an additional boundary condition 
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to the curl-curl equations. They introduced Lagrange multipliers (which axe identical 
to the “dummy” variables used by the authors for other reason, see [18] ) associated 
to the divergence constrains and reformulated the full Maxwell equations as a con- 
strained variational problem. Then they applied the well-developed mixed Galerkin 
finite element method in fluid mechanics with the use of non-equal-order elements 
to solve the problem. This method was one of mathematically solid approaches to 
deal with the time-domain Maxwell equations. 

Boyse and Seidl [5] also realized that node-based methods based on the double 
curl equations suffers from spurious modes. They pointed out that such a spurious 
solution is accompanied by ill-conditioned finite element matrices, non-convergent 
iterative methods, and gross errors in the field calculation. They proposed a hy- 
brid finite element method (HFEM) which combines the node-based FEM and the 
method of moments (MOM). They used a scalar and vector potential formulation 
which eliminates spurious solutions. Moreover, the method based on the Helmholtz 
equation, such as the one used by Bristeau et al. [6] and Dabaghi [11], is also free 
of spurious modes. 

The edge element method [26] is becoming popular for the time domain calculation. 
The author’s opinion regarding edge elements is expressed in [18], and thus will not 
be repeated here. 

Recently, Jiang et al. [15, 18] pointed out that: 


1. the principal part, which consists of the first-derivative terms, of the full 
Maxwell equations in their static, time-harmonic and time-discretized forms, 
contains two div-curl systems; 

2. by introducing a “dummy” variable (this treatment was first proposed by 
Chang and Gunzburger [9] ) which is identically zero, it can be shown that the 
div-curl system is elliptic and properly determined, i.e., it is not “overdeter- 
mined” . Consequently the full Maxwell equations are not “overdetermined” , 
and the divergence equations are not redundant; 

3. the inclusion of the divergence-free constrains is necessary — to guarantee the 
uniqueness of the solution in stationary cases, to improve the accuracy of the 
numerical solution in time- varying cases, and to exclude an infinite degenerate 
eigenvalue in time-harmonic cases; 
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4. the satisfaction of the divergence equations can be made easy by using the 
least-squares method. 


The objectives of this paper are (1) to provide an implicit node- based LSFEM for 
numerical solution of time-dependent full Maxwell equations; (2) to demonstrate 
that the divergence equations should not be disregarded for time-domain problems; 
and (3) to supplement our previous paper [18] with time-dependent examples. 

The LSFEM is based on the min imiz ation of the £2 norm of the residuals of first- 
order systems. As the name implies, the spatial discretization is achieved by using 
the finite element method. The LSFEM is a universal method for numerical solution 
of all types of partial differential equations. Further information on the LSFEM can 
be found in Bochev and Gunzburger [4], Fix and Rose [13], Jiang et al. [14], Jiang 
and Povinelli [16, 17], Krizek and Neittaanmaki [20] and Lager and Mur [22] among 
others. 

The necessity of implicit schemes for the approximate solution of Maxwell equations 
have been clearly explained by Adam et al. [1] and Monk [24]. The attraction of 
imp1ir.it time marching schemes is the avoidance of unnecessarily small time steps 
required by explicit schemes due to the presence of a few small elements in the mesh. 
The implicit scheme presented in this paper is also inexpensive, since the LSFEM 
always leads to a symmetric positive definite system of algebraic equations which 
can be effectively solved by a matrix-free conjugate gradient method. For scattering 
problems, only a few iterations are needed to advance one time step (more details 
are provided in Section 3). 

In the next section the governing equations and the associated boundary conditions 
will be presented in detail. The Maxwell equations are naturally of first order and 
are thus particularly suitable for the LSFEM. The governing equations are first dis- 
cretized in time using the Crank-Nicolson scheme, then the LSFEM is applied to 
obtain numerical solutions. The divergence-free conditions are incorporated in the 
least-squares functional. The discretization procedure is elaborated in Section 3. 
It will become apparent that the present scheme does not require any special tech- 
niques, such as the use of a staggered grid, non-equal-order elements and upwinding, 
thus, it can be implemented into any existing finite element code in a straight for- 
ward manner. The description of the test cases and the numerical results are given 
in Section 4. Three examples are shown. The first is the polarization by a perfectly 
conducting circular cylinder at various incident wavelengthes for both transverse- 
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electric (TE) and transverse-magnetic (TM) waves. The other tests considered axe 
the scattering by a square cylinder and a NACA0012 airfoil. In order to illustrate 
the importance of the inclusion of divergence-free equations, for the circular cylin- 
der problem, the results of the LSFEM based on the full Maxwell equations are 
compared with those based on the solution of only the two curl equations, which 
clearly indicate the occurrence of spurious solutions with large divergence near the 
scatterers when the divergence-free constraints are neglected. The contours of the 
spurious solutions observed here also closely resemble those produced by a FD-TD 
approach [40]. 


2 Governing Equations 

2.1 The Maxwell Equations 

In this paper, the time-dependent Maxwell Equations in free space are of interest, 
which are written as: 


V • E = 0, 

Gauss’ electric law 

(1) 

V • H = 0, 

Gauss’ magnetic law 

(2) 

eE* - V x H = 0, 

Ampere’s law 

( 3 ) 

lxH t + V x E = 0, 

Faraday’s law 

( 4 ) 


where E and H are intensities of the electric and magnetic fields respectively; and the 
physical constants e and /z are the dielectric permittivity and magnetic permeability 
of free space. 

We note that there are eight equations and six unknown variables in Eqs. (1-4). 
It has been shown by Jiang et al. [18] that these systems are not overdetermined 
and that the divergence equations, Eqs. (1) and (2), are not redundant and can be 
satisfied by the least-squares method. 

The objective is to determine the scattered field due to the existence of perfect 
conducting bodies exposed to plane incident waves in free space. For such cases the 
field intensity variables are conveniently split into those of the incident field and of 
the scattered field, such as 

E = E’ + E*, (5) 
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H 1 + IP, 


and H = 


(6) 


where the superscript “i” denotes the incident field; and “s” denotes the scattered 
field. 


By using the above scattered formulation, the incident waves are allowed to prop- 
agate in their analytical form. The scattered waves are determined by solving the 
Maxwell equations with appropriate boundary conditions on the surface of the scat- 
terer and in the fax-field. 


2.2 Boundary Conditions 

On the surface of a perfectly conducting body, the tangential components of the 
electric intensity and the normal component of the magnetic intensity axe both 
zero, i.e., 

n x E = 0 and n-H = 0, (7) 

where n is the outward directed surface normal. 

When a scattered field formulation is used, the above becomes: 


nxE’ = -nxE' and n • IP = — n-H\ (8) 

Since the Maxwell equations are now solved in a finite domain, proper far-field 
boundary conditions have to be applied to ensure that waves will not artificially 
reflect back and contaminate the near-field solution. Here we adopt the following 
absorbing boundary condition: 

1 d<f> d<f> 

c dt + dn 

where n is the outwaxd directed normal of the far-field boundary, (f> is any component 
of E* and IP, and c = l/y/eji is the wave velocity. This absorbing boundary 
condition can be easily implemented in the current version of our general-purpose 
LSFEM code. 

It should be noted that in electromagnetic scattering problems normally only the 
solution in the near-field is of interest. Indeed the radar cross section (RCS) calcu- 
lation requires only the field intensities on the surface of the scatterer (see Section 
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4). Because the absorbing boundary condition given in Eq. (9) is only an approxi- 
mation of the true absorbing boundary condition, the computational domain needs 
to be larger than those using higher-order absorbing conditions. However, at the 
far-field it is only required that the waves propagate out and do not reflect back to 
contaminate the near field solutions, the accuracy of the solution there is of no inter- 
est. For this reason it is most efficient to refine the mesh locally near the surface of 
the scatterer and to use a coarse mesh in the far-field. The use of such non-uniform 
grids is easily accommodated in the present scheme since the calculation is carried 
out totally on an unstructured mesh. Due to the use of a non-uniform grid and large 
elements near the far-field boundary, the elements placed in the part of the com- 
putational domain far from the scatterer surface represent only a small proportion 
of the total element number. Thus the present approach still provides a very good 
appro xim ation of the field variables while maintaining efficiency. We will later show 
in Section 4 that satisfactory results are obtained using such non-uniform meshes. 
We remark that this refinement strategy is not new. Ambrosiano et al. [2] and Vinh 
et al. [40] have used grids with large size variation and obtained satisfactory RCS. 
We also want to emphasize that this concept of local refinement differs from that 
based on a uniform distribution of error, which is widely used in solid mechanics 
and other fields. 


2.3 Transverse-Magnetic (TM) and Transverse-Electric (TE) 
Forms in Two Dimensions 


If the electromagnetic fields are two-dimensional, all derivatives with respect to the 
third dim ension become zero, and some components of the field intensity variables 
vanish. There exist two sets of distinct modes, that is, the Transverse- Magnetic 
(TM) mode and Transverse-Electric (TE) mode. 

For TM waves, we have 

H z = 0, £* = 0, E y = 0. (10) 


The governing equations, Eqs. (1-4), become: 


dE z dH v dH x 


V- 


dt 

dH, 


+ 


dx 

dE z 


dt dy 


+ dy 
= 0 , 


= 0 , 


( 11 ) 

(12) 
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dHy 6E Z n 

“at dx = °’ 

(13) 


SH, SH , 

(14) 


x | y n 

dx dy ~ ' 

Whereas for TE waves, 

E z =Q,H X = 0, H y = 0, 

(15) 

and the governing equations 

are simplified as: 



dH z dE y dE z n 

(16) 


“ at + ax ~ a y = 0 ' 


dE . SH, „ 
c a, dy =0 ' 

(17) 


e dE v + 8H. _ Q 
dt dx 

(18) 


dE x dE y 
dx dy 

(19) 


3 Discretization of the Governing Equations 


Equations (1-4) axe first discretized in time, using the Crank-Nicolson scheme. This 
scheme provides second order accuracy in time and is known to be nondissipative 
when used in conjunction with central space differences. Although this scheme is 
unconditionally stable and thus poses no limitation on the time step, the dispersive 
error becomes large if very large time step is used. Thus, care should be taken when 
choosing the size of the time step. To maintain the accuracy in the time domain, 
we typically use a time step corresponding to a CFL number between 1 and 2 at 
the smallest element. 


After applying time discretization to Eqs. (1-4), we have, 


A»- 


E n+1 _ E n 

: At 

H «+i _ H" 

At 


V • E n+1 = 0 

(20) 

V ■ I T +1 = 0 

(21) 

- x (H n + H n+1 ) = 0 

(22) 

+ -V x (E n + E n+1 ) = 0 
£ 

(23) 
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The two- dim ensional TM and TE cases axe presented in the same form, with the 
corresponding zero terms dropped. 

Eq. (24) is then discretized in space following the standard LSFEM procedure [16]. 
This results in a set of algebraic equations which can be written as: 

KU = F (26) 
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In the above, it is implied that the following standard finite element approximation 
to the unknown variables u has already been introduced, i.e., 


uku = NU (27) 

where u is the finite element approximation to u and U are the nodal values; N are 
the finite element shape functions. 

The matrix K in Eq. (26) is defined by the matrices A and A; in Eq. (24) and 
shape functions N in Eq. (27) as: 

K = / [L(N)] r L(N)<fl2 (28) 

J n 

in which 

L(N) = AN + AiN rf 

i=l,3 

and [L(N)] t denotes the transpose of the matrix L(N); 12 is the domain of the 
problem. 

The right hand side vector in Eq. (26) is defined by: 

F = / [L(N)] T fdft (29) 

J n 

where f is the right hand side of equation (24). 

At this stage it becomes clear that once the governing equation set is written in the 
form of the standard first-order system, i.e., Eq. (24), the LSFEM can be easily 
incorporated into any existing finite element code in a straight forward manner. 
One advantage of using the LSFEM is that the computer code can be constructed 
in a general-purpose setting for the system of equations in Eq. (24) in such a way 
that the only required programming for application to a new problem is to write a 
subroutine to supply the matrices A, A j and f, therefore the tedious work is reduced 
to a minimum. 

We have used a Jacobi-preconditioned conjugate gradient method to solve the re- 
sulting algebraic equation system, Eq. (26). By using this method, the iterations 
ran be performed in such a matrix-free manner that there is no need to assemble 
the global matrix K or to form the element matrices [8]. This has greatly reduced 
the requirement for computer memory. We have found that five conjugate gradient 
iterations are enough to guarantee the accuracy of the solution in our calculations. 
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The use of such a small number of iterations can be attributed to the fact that 
when the time step is small (corresponding to a CFL number between 1 and 2), the 
solution in the previous time step serves as a very good initial guess for the next 
step, thus the conjugate gradient method needs fewer iterations than the case when 
a larger time step is used. Therefore, although the current method takes the form of 
an implicit method, its computational speed is competitive with an explicit method. 


4 Numerical Examples 


Several test cases of the scattering of a plane incident wave by perfect conduct- 
ing bodies will be presented. Since the purpose of this paper is to validate the 
present scheme and show the importance of the divergence equations, the numerical 
examples will be confined to two dimensional problems. We believe these cases suf- 
ficiently demonstrate our view points. We have successfully carried out preliminary 
calculations for three dimensional problems and have recently reported these results 

[43]. 

The incident sinusoidal wave is given by: 


E z ' x — i£osinfc(xcosa 4- ysina — ct), 


H 


*.« _ 


H y,i = 


E z,t sin a 

fic 

E z,l cosa 


fic 


for TM case; and 


H z,t = Ho sin k(x cos a + y sin a — ct ), 


E x,i 


H z,t sin a 

£C 


E v,i 


H z ’ 1 cos a 


ec 

for TE case. Here c is the wave velocity defined in Section 2; k is the wave 
a is the angle of incidence. 


number; 


The numerical results will be presented in terms of the surface current (for TM 
case only) and radar cross section (RCS). The calculation of these requires the field 
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variables in the frequency domain, thus a Fourier transform is performed to transfer 
the field variables from the time domain to the frequency domain. 

The surface current J s for TM case is defined as 

jn x H,| 

' |H(| ’ 

where H* and H t are the incident and total magnetic fields in the frequency domain, 
respectively. 

The radar cross section (RCS) for TM wave, S TM , is defined as: 

< 30 > 

where E* z and E\ are the scattered and incident electric fields, respectively, in the 
frequency domain; p is the distance from the center of the scattering object; and 
6 is the angle of observation. Note that the calculation of S TM from the above 
expression requires knowledge of E* at an infinite distance from the scatterer. This 
is not available since the Maxwell equations are solved in a finite domain. However, 
the far-field solution can be obtained from the distribution on a near-field surface 
(the most convenient one being the surface of the scatterer) by a transformation 
using Green’s function, as described by Shankar et al. [33]. This transformation is 
adopted in this paper to obtain the RCS. 

For TE wave, the RCS S TB is defined as: 

S„(#)=Hm2xp|^^| 2 . (31) 

The calculation of S TE follows the same procedure as TM cases. 

Bilinear quadrilateral elements for all variables are used in all cases. The meshes 
are generated using the unstructured mesh generator developed by Zhu et al. [45], 
based on the advancing front technique of Peraire et al. [29]. 

4.1 Circular Cylinder 

The first example deals with a scattering field due to a circular perfect conducting 
circular cylinder. Both TM and TE cases at various wavelengthes are considered. 
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The definition of the problem is given in Fig. 1. This problem allows an analytical 
solution in a series form, see, e.g., Umashankar and Taflove [39]. Due to the sym- 
metry of this problem, only a half of the domain is used. The computational mesh 
consists of 2538 nodes and 2435 elements, with 142 divisions on the cylinder surface 
(Fig. 2). Figs. 3-4 show the calculated RCS. Fig. 5 shows the calculated surface 
current for TM cases and compares the LSFEM results with the analytical solutions 
which are calculated from the series form presented in [39] using Mathematica [41]. 
The instantaneous distributions of the scattered field variables are shown in Figs. 
6-7. 


4.2 Square Cylinder 

The second example concerns scattering by a square cylinder defined in Fig. 8. 
The TM case with ka = 1 is considered. The computational mesh consists of 2749 
bilinear quadrilateral elements with 2874 nodes (Fig. 9). Fig. 10 shows the surface 
current distribution. Comparison is made with the results by Shankar et al. [34] 
and that obtained by Umashankar and Taflove [39] using the method of moments 
(MOM). The instantaneous E z contours axe shown in Fig. 11. 


4.3 NACA0012 Airfoil 


Preliminary results are presented here for two cases of scattered TM wave due to a 
NACA0012 airfoil (Fig. 12). The angles of incidence are 0 and 90 degrees; and the 
values of ka are 10 and 107r, respectively. Fig. 13 shows the computational mesh, 
which consists of 1782 bilinear elements and 1912 nodes. There are 156 divisions on 
the surface of the airfoil. Fig. 14 shows the contours of Ez for the 0 degree incidence 
case. Figs. 15 and 16 show the computed RCS. The case of 90 degree incidence is 
compared with those by Shankar et al. [34] and Vinh et al. [40]. Generally, fair 
agreement is observed. Note that we have used a very coarse mesh and a very small 
domain. 
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4.4 Importance of the Divergence Equations 

To illustrate the influence of the inclusion of divergence equations, we present here 
a comparison of two groups of numerical solutions, both obtained by using the 
LSFEM: the first group is obtained by solving the full Maxwell equations, while the 
second group is obtained by solving solely the two curl equations. Here the case of 
TM waves scattered by a circular cylinder is considered. The values of ka used here 
are ka = 1 and ka = 10. Fig. 17 shows a comparison of the contours of the field 
variables for ka = 1 between the group 1 and the group 2. It is clearly revealed 
that without enforcing the divergence condition significant spurious patterns exist 
for H x and H y near the surface. Further calculation shows that indeed the group 2 
solution possesses very large divergence of H in the vicinity of the scatterer surface. 
The absolute value of V • H reaches its maximum of 21.28 in an element adjacent 
to the surface. On the other hand, the maximum value of V • H for the group 1 is 
0.054, which occurs near the outer boundary. The divergence-free condition for E 
is automatically satisfied, therefore almost identical E z distributions are obtained 
with the two formulations. Fig. 18 shows the field variables for ka = 10. Again 
a non-physical solution is found for the group 2 solution near the surface. These 
non-physical patterns are similar to those in Figs. 2-b and 2-c of [40]. 

We should remark here that although group 1 and group 2 solutions differ signifi- 
cantly, very little difference was observed in the RCS calculations based on the two 
solutions. This is due to the fact that the erroneous part of the solution is static, 
and the static field does not radiate. Definitely, the erroneous solutions of the field 
intensities caused by the neglect of the divergence conditions cannot be accepted, for 
example, in plasma computations and in the computation of impedance in digital 
circuits and electronic packaging. 


5 Concluding Remarks 


A node-based least-squares finite element method has been presented for the time- 
dependent Maxwell equations for scattering problems. Divergence-free conditions 
have been included in the least-squares functional thus eliminating the existence 
of spurious modes. This technique is also expected to be useful for the solution of 
Vlasov-Maxwell equations in plasma physics and the simulation of magnetohydrody- 
namics (MHD), and for the modeling of microwave integrated circuit components. 
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The accuracy of the scheme has been verified by various benchmark tests. The 
computation is carried out on a totally unstructured mesh which has a rather large 
distribution of grid sizes. This provides a great amount of flexibility in dealing with 
complex geometries. The scheme is implicit and unconditionally stable, therefore, 
it offers significant advantages in handling electromagnetic problems whose element 
size varies several orders of magnitude across the problem domain. In those situa- 
tions the time step will no longer be determined by the smallest size as in explicit 
schemes. Although this method is implicit, due to the use of the matrix-free con- 
jugate gradient method, there is no matrix inversion, and thus the computational 
cost and memory requirements are competitive with explicit approaches. 

In this report, we have demonstrated that the omission of the divergence equations in 
the numerical procedure produces non-physical solutions for scattered field problems, 
we have therefore shown that the common belief that the divergence equations are 
redundant is a serious misconception. The possibility of producing spurious solutions 
from time-domain methods should not be overlooked, even though it does not affect 
the computed RCS. 
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Figure 1 Scattering by a circular cylinder: problem definition 



Figure 2 Scattering by a circular cylinder: finite element mesh 





















Figure 6 Scattered fields for the circular cylinder 

( E z , H x , H y contours and H vectors, TM mode, ka = 5) 
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Figure 7 Scattered fields for the circular cylinder 

( E z , H x , H y contours and H vectors, TM mode, ka = 1) 
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Figure 8 Scattering by a square cylinder: problem definition 



Figure 9 Scattering by a square cylinder: finite element mesh 
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Figure 11 Scattered E z contours for the square cylinder 
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Figure 12 Scattering by a NACA0012 aerofoil: problem definition 



Figure 13 Scattering by a NACA0012 aerofoil: finite element mesh 



Figure 14 Scattered E z contours for the NACA0012 aerofoil 
(0 degree incidence, ka = 10 ) 
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NACA0012, 0 degree incidence 



Figure 15 RCS for the NACA0012 airfoil (0 degree incidence, ka = 10 ) 



Figure 16 RCS for the NACA0012 airfoil (90 degree incidence, ka - 10tt ) 
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